Volume 2: The Fourier Transform.

<Name>
<Class>
<Date>

Part 1: The Discrete Fourier Transform

In [1]:
from matplotlib import pyplot as plt
from scipy.io import wavfile
from scipy import signal
import numpy as np
import IPython
from scipy import fftpack
import math
import imageio
In [2]:
plt.rcParams["figure.dpi"] = 300             # Fix plot quality.
plt.rcParams["figure.figsize"] = (12,3)      # Change plot size / aspect (you may adjust this).
In [3]:
class SoundWave(object):
    """A class for working with digital audio signals."""

    # Problem 1.1
    def __init__(self, rate, samples):
        """Set the SoundWave class attributes.

        Parameters:
            rate (int): The sample rate of the sound.
            samples ((n,) ndarray): NumPy array of samples.
        """
        self.rate = rate
        self.samples = samples

    # Problems 1.1 and 1.7
    def plot(self, dft = False):
        """Plot the graph of the sound wave (time versus amplitude)."""
        
        #plot regardless of input
        ax1 = plt.subplot(211)
        domain = np.linspace(0, len(self.samples)/self.rate, len(self.samples))
        
        ax1.plot(domain, self.samples)
        ax1.set_ylim(-32768, 32767)
        plt.xlabel("Seconds")
        plt.ylabel("Magnitude")
        if not dft:
            plt.show()
            
        #Only plot if dft is True
        if dft:
            #calculate
            c_k = abs(fftpack.fft(self.samples))
            n = len(c_k)
            k_indices = np.arange(n)
            frequency = k_indices * self.rate/n
            #plot
            ax2 = plt.subplot(212)
            ax2.plot(frequency, c_k)
            ax2.set_xlim((0,self.rate/2))
            plt.show()
            

    # Problem 1.2
    def export(self, filename, force=False):
        """Generate a wav file from the sample rate and samples. 
        If the array of samples is not of type np.int16, scale it before exporting.

        Parameters:
            filename (str): The name of the wav file to export the sound to.
        """
        if force or self.samples.dtype != np.int16:
            #Scaling
            scaled_samples = np.int16((self.samples * 32767.0)/max(abs(self.samples)))
            wavfile.write(filename, self.rate, scaled_samples)
        else:
            wavfile.write(filename, self.rate, self.samples)
    
    # Problem 1.4
    def __add__(self, other):
        """Combine the samples from two SoundWave objects.

        Parameters:
            other (SoundWave): An object containing the samples to add
                to the samples contained in this object.
        
        Returns:
            (SoundWave): A new SoundWave instance with the combined samples.

        Raises:
            ValueError: if the two sample arrays are not the same length.
        """
        if len(self.samples) != len(other.samples):
            raise ValueError("sample lengths not the same size")
        else:
            #add matricies elementwise
            new_samples = self.samples + other.samples
            return SoundWave(self.rate, new_samples)

    # Problem 1.4
    def __rshift__(self, other):
        """Concatentate the samples from two SoundWave objects.

        Parameters:
            other (SoundWave): An object containing the samples to concatenate
                to the samples contained in this object.

        Raises:
            ValueError: if the two sample rates are not equal.
        """
        if self.rate != other.rate:
            raise ValueError("rates are not equal")
        else:
            #concatenate matricies
            new_samples = np.concatenate((self.samples,other.samples))
            return SoundWave(self.rate, new_samples)
    
    # Problem 2.1
    def __mul__(self, other):
        """Convolve the samples from two SoundWave objects using circular convolution.
        
        Parameters:
            other (SoundWave): An object containing the samples to convolve
                with the samples contained in this object.
        
        Returns:
            (SoundWave): A new SoundWave instance with the convolved samples.

        Raises:
            ValueError: if the two sample rates are not equal.
        """
        if self.rate != other.rate:
            raise ValueError("Rates not equal")
        else:
            if len(self.samples) < len(other.samples):
                self.samples = np.pad(self.samples,(0,len(other.samples) - len(self.samples)),'constant')
            elif len(self.samples) > len(other.samples):
                other.samples = np.pad(other.samples,(0,len(self.samples) - len(other.samples)),'constant')
            #convolve arrays
            convolution = fftpack.ifft(fftpack.fft(self.samples) * fftpack.fft(other.samples))
            return SoundWave(self.rate, convolution)

    # Problem 2.2
    def __pow__(self, other):
        """Convolve the samples from two SoundWave objects using linear convolution.
        
        Parameters:
            other (SoundWave): An object containing the samples to convolve
                with the samples contained in this object.
        
        Returns:
            (SoundWave): A new SoundWave instance with the convolved samples.

        Raises:
            ValueError: if the two sample rates are not equal.
        """
        if self.rate != other.rate:
            raise ValueError("Rates not equal")
        else:
            n = len(self.samples)
            m = len(other.samples)
            length = n + m - 1
            a = np.log2(length)
            a = math.ceil(a)
            
            self_sam = np.pad(self.samples,(0,2**a - n),'constant')
            other_sam = np.pad(other.samples,(0,2**a - m),'constant')
            
            #convolve arrays
            convolution = fftpack.ifft(fftpack.fft(self_sam) * fftpack.fft(other_sam))
            #take first n+m-1 elements
            convolution = convolution[:length]
            
            return SoundWave(self.rate, convolution)

    # Problem 2.4
    def clean(self, low_freq, high_freq):
        """Remove a range of frequencies from the samples using the DFT. 

        Parameters:
            low_freq (float): Lower bound of the frequency range to zero out.
            high_freq (float): Higher boound of the frequency range to zero out.
        """
        n = len(self.samples)
        k_low = math.floor(low_freq * n / self.rate)
        k_high = math.ceil(high_freq * n / self.rate)
        
        dft = fftpack.fft(self.samples)
        dft[k_low:k_high] = 0 
        dft[n - k_high: n - k_low] = 0
        
        self.samples = fftpack.ifft(dft).real

Problem 1.1

  • Implement SoundWave.__init__().
  • Implement SoundWave.plot().
  • Use the scipy.io.wavefile.read() and the SoundWave class to plot tada.wav.
In [4]:
rate, samples = wavfile.read("tada.wav") #pull info from file
a = SoundWave(rate,samples)
a.plot()

Problem 1.2

  • Implement SoundWave.export().
  • Use the export() method to create two new files containing the same sound as tada.wav: one without scaling, and one with scaling (use force=True).
  • Use IPython.display.Audio() to embed the original and new versions of tada.wav in the notebook.
In [5]:
IPython.display.Audio(filename = 'tada.wav') #read file for comparison
Out[5]:
In [6]:
rate, samples = wavfile.read("tada.wav")
tada_ = SoundWave(rate,samples)
tada_.export('tada_not_scaled.wav') #export file, should be same as previous
IPython.display.Audio(filename = 'tada_not_scaled.wav')
Out[6]:
In [7]:
tada_.export('tada_scaled.wav', True) #export scaled file, should be louder than previous
IPython.display.Audio(filename = 'tada_scaled.wav')
Out[7]:

Problem 1.3

  • Implement generate_note().
  • Use generate_note() to create an A tone that lasts for two seconds. Embed it in the notebook.
In [8]:
def generate_note(frequency, duration):
    """Generate an instance of the SoundWave class corresponding to 
    the desired soundwave. Uses sample rate of 44100 Hz.
    
    Parameters:
        frequency (float): The frequency of the desired sound.
        duration (float): The length of the desired sound in seconds.
    
    Returns:
        sound (SoundWave): An instance of the SoundWave class.
    """
    Rate = 44100 #code from textbook
    num_samples = duration * Rate
    domain = np.linspace(0,duration, num_samples)
    samples = np.sin(2*np.pi*domain*frequency)
    return SoundWave(Rate, samples)
In [9]:
a = generate_note(440, 2) #A note

IPython.display.Audio(rate=a.rate, data=a.samples)
Out[9]:

Problem 1.4

  • Implement SoundWave.__add__().
  • Generate a three-second A minor chord (A, C, and E) and embed it in the notebook.
  • Implement SoundWave.__rshift__().
  • Generate the arpeggio A$\,\rightarrow\,$C$\,\rightarrow\,$E, where each tone lasts one second, and embed it in the notebook.
In [10]:
a = generate_note(440, 3)
c = generate_note(523.25, 3)
e = generate_note(659.25,3)
cord = a+c+e #notes added elementwise
IPython.display.Audio(rate=cord.rate, data=cord.samples)
Out[10]:
In [11]:
a = generate_note(440, 1)
c = generate_note(523.25, 1)
e = generate_note(659.25,1)
arpeggio = a>>c>>e #notes concatenated
IPython.display.Audio(rate=arpeggio.rate, data=arpeggio.samples)
Out[11]:

Problem 1.5

  • Implement simple_dft() with the formula for $c_k$ given below.
  • Use np.allclose() to check that simple_dft() and scipy.fftpack.fft() give the same result (after scaling).
$$ c_k = \frac{1}{n}\sum_{\ell=0}^{n-1} f_\ell e^{-2 \pi i k \ell\, /\, n} $$
In [12]:
def simple_dft(samples):
    """Compute the DFT of an array of samples.

    Parameters:
        samples ((n,) ndarray): an array of samples.
    
    Returns:
        ((n,) ndarray): The DFT of the given array.
    """
    n = len(samples) #code from textbook
    m = np.arange(n).reshape(n,1)
    W = np.exp((-2j*np.pi/n)*m@m.T)
    return W@samples/n
In [13]:
"""
n = 8192 #power of 2
samples = np.random.random(n)
simple = simple_dft(samples)
sci = fftpack.fft(samples)
#compare textbook dft with scipy dft
np.allclose(n*simple,sci)
"""
Out[13]:
'\nn = 8192 #power of 2\nsamples = np.random.random(n)\nsimple = simple_dft(samples)\nsci = fftpack.fft(samples)\n#compare textbook dft with scipy dft\nnp.allclose(n*simple,sci)\n'

Problem 1.6

  • Implement simple_fft().
  • Generate an array of $8192$ random samples and take its DFT using simple_dft(), simple_fft(), and scipy.fftpack.fft(). Print the runtimes of each computation.
  • Use np.allclose() to check that simple_fft() and scipy.fftpack.fft() give the same result (after scaling).
In [14]:
def simple_fft(samples, threshold=1):
    """Compute the DFT using the FFT algorithm.
    
    Parameters:
        samples ((n,) ndarray): an array of samples.
        threshold (int): when a subarray of samples has fewer
            elements than this integer, use simple_dft() to
            compute the DFT of that subarray.
    
    Returns:
        ((n,) ndarray): The DFT of the given array.
    """
    #code from textbook
    n = len(samples) 
    if n <= threshold:
        return simple_dft(samples)
    else:
        even = simple_fft(samples[::2])
        odd = simple_fft(samples[1::2])
        w = np.exp(-2j*np.pi/n * np.arange(n))
        m = n//2
        first_sum = even + w[:m] * odd
        second_sum = even + w[m:] * odd
        return 0.5 *np.concatenate([first_sum, second_sum])
In [15]:
'''
n = 8192
samples = np.random.random(n)
#time computation times
%time n*simple_dft(samples)
%time n*simple_fft(samples)
%time fftpack.fft(samples)
'''
Out[15]:
'\nn = 8192\nsamples = np.random.random(n)\n#time computation times\n%time n*simple_dft(samples)\n%time n*simple_fft(samples)\n%time fftpack.fft(samples)\n'

Problem 1.7

  • Modify SoundWave.plot() so that it accepts a boolean. When the boolean is True, take the DFT of the stored samples and plot (in a new subplot) the frequencies present on the $x$-axis and the magnituds of those frequences on the $y$-axis. Only the display the first half of the plot, and adjust the $x$-axis so that it correctly shows the frequencies in Hertz.
  • Display the plot of the DFT of an A tone.
  • Display the plot of the DFT of an A minor chord.
In [16]:
a.plot(True)
In [17]:
cord.plot(True)

Problem 1.8

Use the DFT to determine the individual notes that are present in mystery_chord.wav.

In [18]:
rate, samples = wavfile.read("mystery_chord.wav")
mystery = SoundWave(rate,samples)

IPython.display.Audio(filename = 'mystery_chord.wav')
mystery.plot(True)

# Find the Magnitudes

n = len(mystery.samples)
k_indices = np.arange(n)
freq = (mystery.rate * k_indices) / n
magnitude = np.abs(fftpack.fft(mystery.samples))

# Find the frequencies
sort = np.argsort(magnitude[:n//2])[::-1]
note_1 = freq[sort[0]]
note_2 = freq[sort[1]]
note_3 = freq[sort[2]]
note_4 = freq[sort[3]]

print(note_1, 'A')
print(note_2, 'G')
print(note_3, 'C')
print(note_4, 'D')

440.0 A
784.0 G
523.25 C
587.5 D
In [ ]:

The notes are...

Part 2: Convolution and Filtering.

Problem 2.1

  • Implement SoundWave.__mul__() for circular convolution.
  • Generate 2 seconds of white noise at the same sample rate as tada.wav.
  • Compute the circular convolution of tada.wav and the white noise. Embed the result in the notebook.
  • Append the circular convolution to itself and embed the result in the notebook.
In [19]:
rate_wn = 22050        # Create 2 seconds of white noise at a given rate.
white_noise = np.random.randint(-32767, 32767, rate_wn*2, dtype=np.int16)
noise_sound = SoundWave(rate_wn, white_noise)
IPython.display.Audio(rate=noise_sound.rate, data=noise_sound.samples)
Out[19]:
In [20]:
#convolution
convo = tada_*noise_sound
#append convolution to itself
new_convo = convo >> convo
IPython.display.Audio(rate=new_convo.rate, data=new_convo.samples)

/Users/garretcarver/anaconda3/lib/python3.6/site-packages/IPython/lib/display.py:124: ComplexWarning: Casting complex values to real discards the imaginary part
  data = np.array(data, dtype=float)
Out[20]:

Problem 2.2

  • Implement SoundWave.__pow__() for linear convolution.
  • Time the linear convolution of CGC.wav and GCG.wav using SoundWave.__pow__() and scipy.signal.fftconvolve().
  • Embed the two original sounds and their convolutions in the notebook. Check that the convolutions with SoundWave.__pow__() and scipy.signal.fftconvolve() sound the same.
In [21]:
rate, samples = wavfile.read("CGC.wav")
cgc_sound = SoundWave(rate,samples)

IPython.display.Audio(rate=cgc_sound.rate, data=cgc_sound.samples)
Out[21]:
In [22]:
rate, samples = wavfile.read("GCG.wav")
gcg_sound = SoundWave(rate,samples)
IPython.display.Audio(rate=gcg_sound.rate, data=gcg_sound.samples)
Out[22]:
In [23]:
new_tone = cgc_sound ** gcg_sound
IPython.display.Audio(rate=new_tone.rate, data=new_tone.samples)

/Users/garretcarver/anaconda3/lib/python3.6/site-packages/IPython/lib/display.py:124: ComplexWarning: Casting complex values to real discards the imaginary part
  data = np.array(data, dtype=float)
Out[23]:
In [24]:
scipy_tone = signal.fftconvolve(cgc_sound, gcg_sound)
#IPython.display.Audio(rate=new_tone.rate, data=new_tone.samples)

%time cgc_sound ** gcg_sound
%time signal.fftconvolve(cgc_sound, gcg_sound)

CPU times: user 173 ms, sys: 20 ms, total: 193 ms
Wall time: 104 ms
CPU times: user 62.8 ms, sys: 4.83 ms, total: 67.6 ms
Wall time: 38.2 ms
Out[24]:
<__main__.SoundWave at 0x1c2164a860>

Problem 2.3

Use SoundWave.__pow__() or scipy.signal.fftconvolve() to compute the linear convolution of chopin.wav and balloon.wav. Embed the two original sounds and their convolution in the notebook.

In [25]:
rate, samples = wavfile.read("chopin.wav")
choping = SoundWave(rate,samples)
IPython.display.Audio(rate=choping.rate, data=choping.samples)
Out[25]:
In [26]:
rate, samples = wavfile.read("balloon.wav")
balloon = SoundWave(rate,samples)
IPython.display.Audio(rate=balloon.rate, data=balloon.samples)
Out[26]:
In [27]:
choping_balloon = choping ** balloon
IPython.display.Audio(rate=choping_balloon.rate, data=choping_balloon.samples)

/Users/garretcarver/anaconda3/lib/python3.6/site-packages/IPython/lib/display.py:124: ComplexWarning: Casting complex values to real discards the imaginary part
  data = np.array(data, dtype=float)
Out[27]:

Problem 2.4

  • Implement SoundWave.clean().
  • Clean noisy1.wav by filtering out frequencies from $1250$-$2600$ Hz. Embed the original and the cleaned versions in the notebook.
  • Clean noisy2.wav. Embed the original and the cleaned versions in the notebook.
In [28]:
rate, samples = wavfile.read("noisy1.wav")
noisy1 = SoundWave(rate,samples)

noisy1.clean(1250, 2600)

IPython.display.Audio(rate=noisy1.rate, data=noisy1.samples)
Out[28]:
In [29]:
rate, samples = wavfile.read("noisy2.wav")
noisy2 = SoundWave(rate,samples)

noisy2.clean(1250, 4500)

IPython.display.Audio(rate=noisy2.rate, data=noisy2.samples)
Out[29]:
In [ ]:

Problem 2.5

  • Clean vuvuzela.wav by filtering bad frequencies out of the left and right channels individually.
  • Recombine the left and right channels and embed the result in the notebook.
In [30]:
rate, samples = wavfile.read("vuvuzela.wav")
vuv1 = SoundWave(rate,samples[:,0])
vuv2 = SoundWave(rate,samples[:,1])

vuv1.clean(200,500)
vuv2.clean(200,500)
vuv = SoundWave(rate, np.vstack((vuv1.samples, vuv2.samples)))

IPython.display.Audio(rate=vuv.rate, data=vuv.samples)
Out[30]:
In [ ]:
In [ ]:
In [ ]:

Problem 2.6

  • Clean up license_plate.png so that the year printed on the sticker in the bottom right corner of the plate is legible.
  • Display the original and cleaned images.
In [31]:
image = imageio.imread("license_plate.png")
plt.imshow(image, cmap = 'gray')
plt.show()
In [32]:
im_dft = fftpack.fft2(image)
plt.imshow(np.log(np.abs(im_dft)), cmap="gray")
plt.show()

/Users/garretcarver/anaconda3/lib/python3.6/site-packages/scipy/fftpack/basic.py:160: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  z[index] = x

The year on the sticker is...2013

In [35]:
mean = np.mean(im_dft)
im_dft[25:45, 90:110] = 350
im_dft[178:184, 328:336] = 330
im_dft[145:155, 225:240] = 320
im_dft[150:250, 300:450] = 350

plt.imshow(np.log(np.abs(im_dft)), cmap="gray")

new_samples = fftpack.ifft2(im_dft).real
plt.imshow(new_samples, cmap = 'gray')
plt.show()

/Users/garretcarver/anaconda3/lib/python3.6/site-packages/scipy/fftpack/basic.py:160: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  z[index] = x